
% Here we want to find Q. (June 11th)

function Q = findQ(tempS, I, p, Nmax)
p1 = p(1); p2 = p(2); p3 = p(3);
Q = zeros(I,Nmax+1);     % Notice: both Q and R's dimention is I * (Nmax+1).
for i = 1:8
    Pre_i = sum(tempS(i,:));
    for m = 1 : (Pre_i+1)
        % Then, how to find the combination?
         for type1 = 0:tempS(i,1)
            for type2 = 0:tempS(i,2)
                for type3 = 0:tempS(i,3)
                    if type1+type2+type3 == m-1
                    	Q(i,m) = Q(i,m)+ ...
                            nchoosek(tempS(i,1),type1)*(p1^(type1))*((1-p1)^(tempS(i,1)-type1))*...
                            nchoosek(tempS(i,2),type2)*(p2^(type2))*((1-p2)^(tempS(i,2)-type2))*...
                            nchoosek(tempS(i,3),type3)*(p3^(type3))*((1-p3)^(tempS(i,3)-type3));
                    end
                end
            end
         end
    end
end
